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ABSTRACT 

An adaptive spectral element method has been developed for the efficient solution of time depen- 
dent partial differential equations. Adaptive mesh strategies that include resolution refinement 
and coarsening by three different methods are illustrated on solutions to the one-dimensional 
viscous Burgers equation and the two-dimensional Navier-Stokes equations for driven flow in a 
cavity. Sharp gradients, singularities and regions of poor resolution are resolved optimally as 
they develop in time using error estimators which indicate the choice of refinement to be used. 
The adaptive formulation presents significant increases in efficiency, flexibility and general ca- 
pabilities for high order spectral methods. 


lr This work was supported in part by the National Aeronautics and Space Administration under NASA 
Contract No. NASl-19480 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, Virginia 23665. 
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1. Introduction 


The adaptive formulation of the spectral element method is aimed at increasing the flexibility 
and range of capabilities of high order spectral methods in general. While spectral methods 
provide highly accurate solutions to partial differential equations governing complex physical 
phenomena, they have been limited to idealized research problems due to their lack of geometric 
flexibility [e.g. 1,2]. Further, while they offer exponential convergence for infinitely smooth 
solutions [3], they have not been useful for problems presenting singularities or thin layers. 
The spectral element method [4] was developed to increase the geometrical flexibility of high 
order spectral methods. In this paper, we present an adaptive spectral element method which 
automatically allocates resolution where it is most needed in an optimal fashion. Singularities 
and thin internal and boundary layers are resolved efficiently as they develop in time. Coupled 
with advances in computer power, the development of an adaptive spectral element method 
presents a tremendous opportunity for high accuracy solutions to “real” engineering problems, 
by reducing the needed computer resources, in terms of storage and cpu time, by several orders 
of magnitude. 

Previous work in nonconforming discretizations [5] and error estimators [6] for the spectral 
element method constituted a first step towards the adaptive formulation. Nonconforming 
discretizations, which allow arbitrary element matchup in the grid, account for a substantial 
savings in resolution over the structured conforming grid spectral element method. A further 
gain in efficiency can be achieved by automating the resolution allocation process. For this 
purpose, error estimators were developed to indicate local elemental error values as well as 
quality of resolution, as measured by decay rates in the solution spectrum. Some less theoretical 
yet very practical issues of an adaptive method are investigated here. A summary of adaptive 
mesh strategies is presented. 

Adaptive mesh capabilities rely on a consistent and efficient rule-based system for refining, 
coarsening and reconstructing a mesh. The components of this rule-based system are determined 
and illustrated by solving the viscous Burgers equation in one dimension. The relative merits of 
changing element size, element position, the number of elements and the order of polynomials 
locally as well as globally is explored. These adaptive refinement options are selected based on 
global mesh optimization criteria as well as local elemental error estimators. The criteria provide 
the necessary constraints on the overall mesh refinement process to ensure global efficiency and 
optimization. 

Adaptive mesh capabilities are illustrated by two examples. The solution of the one- 
dimensional viscous Burgers equation illustrates the adaptive procedure in the presence of 
sharp (low viscosity) and weak (high viscosity) gradients. This newly flexible formulation for 
spectral methods in general, will allow physical problems with sharp gradients, such as com- 
pressible flows with shocks, to be treated automatically, irrespective of the physical application. 
A two-dimensional fluid flow simulation illustrates the ability of the adaptive procedure to treat 
singularities. The full incompressible Navier-Stokes equations are solved for a highly accurate 
simulation of laminar flow in a geometry presenting functional singularities. 


2, Formulation 

In this section the components of the adaptive formulation are described. Together, they 
form a rule-based system designed to produce an optimally efficient and accurate solution to 
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any time-dependent solution of physical phenomena governed by partial differential equations. 
The adaptive formulation relies on an efficient and accurate estimate of the numerical error 
incurred by discretization. The error estimators are the subject of a detailed paper [6]. They 
will be reviewed briefly here. 

2.1 Error Estimator 


Single mesh a posteriori error estimators were developed to estimate the error and indicate 
the quality of resolution on each element. They rely on the calculation and extrapolation of the 
spectrum of the solution on each element to estimate the error and to predict convergence as 
well. Since the spectral element method offers several refinement options, namely increasing the 
polynomial order or increasing the number of elements, it is important that the error estimator 
be able to distinguish which of these refinement options is optimal. The decay rate of the 
spectrum offers information on the quality of the solution. A low decay rate is indicative of 
poor resolution or the presence of a singularity. A high decay rate, on the other hand, indicates 
that the solution is well-resolved. This information is used in the refinement process. 

The error estimate t est used is given in terms of the spectrum a n of the numerical solution 
u/j, defined by the elemental spectral discretization: 


4( r ) = E a t p n(r) 


n=0 


(i) 


where P n is the nth order Legendre polynomial, and N is the discretization order on element 
k. r is the local elemental spatial coordinate. The error estimate is calculated as 
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This is an approximation to the C 2 error. The error in the Ti} norm may be calculated anal- 
ogously. The integral in Equation (2) represents the extrapolation of the spectrum to infinity, 
a measure of the truncation error. The function a(n) is a least squares best fit of the last four 
points of the spectrum to an exponential decay: 


a(n) = ce <T . 


( 3 ) 


The decay rate a indicates poor resolution for a < 1 and good resolution for a > 1. In the 
adaptive process, this information is used to refine by increasing the number of elements and 
increasing the order of the polynomial respectively. The value e e3 t on each element is used to 
decide whether to adapt or not. In practice, we find that these error estimators are very robust 
and quite accurate as shown in [6] and in the following examples. 


2.2 Refinement Criteria 


The refinement decision is based on several criteria. The first and most effective is to 
compare elemental error estimators to a globally acceptable level of error, set once for the 
whole run. The elements with errors over the acceptable level are marked for refinement. The 
elemental error estimators are also compared in a relative manner to neighbouring elements in 
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order to determine which element has the greatest need for increased resolution. This criterion is 
implemented for situations where resolution is limited and only certain elements may be refined. 
It is also used in the no-cost “refinement” option which simply changes element sizes without 
incurring any increase in resolution. Limits and tolerances are imposed at every step of the 
decision process. At each refinement step the minimum and maximum error are calculated. If 
they are both below the acceptable level, no refinement is needed and the calculation proceeds. 
When comparing error estimates between elements, there must be a substantial difference in 
error, usually a factor of two or more, to be able to distinguish which element needs more 
resolution. Further, if any element has a substantially higher error than the minimum error, 
usually a factor of five or more, then it is also marked for refinement. There are also situations 
where one cannot refine. For these we implement limits that prevent refinement from occurring. 
If maximum values of N, the order of the polynomial, K, the number of elements, or NTOTAL, 
the total number of degrees of freedom, are reached refinement is prohibited. Similarly, if 
minimum values of element size or time step size, which decreases with resolution in order to 
satisfy stability conditions, are reached refinement is prohibited. 

The decision to refine by adding elements or increasing the order of the polynomial is 
straightforward. The decision to move elements or to coarsen is not as easy. Coarsening 
implies that the solution does not need as much resolution as it has. While the spectrum decay 
predicts convergence if one adds resolution, it does not have the ability to predict convergence 
if one removes resolution. Theoretically, one might be able to look at the level of the spectrum 
coefficients and estimate what error would be incurred by removing the last coefficient. However, 
this is not very robust, particularly in transient problems. For these reasons, coarsening is 
limited to moving elements. The criterion for moving elements again relies on comparing 
neighbouring elemental error estimates. However, it is difficult to write a democratic algorithm 
which does not destroy the mesh, by having too many elements changing at once. For this 
purpose, a voting system was implemented to mark elements with the grestest need for increased 
“resolution”. The algorithm examines each element and its neighbours, compares their error 
estimates and accordingly votes for the element with the greatest error. After passing through 
the whole grid, a list of the elements with the greatest number of votes is obtained. At the same 
time, an indication of the relative error is saved to serve as a value for shrinking each marked 
element and to indicate from which neighbour the resolution may be taken. Again, a limit is 
imposed to prevent drastic size reductions at the expense of neighbouring elements. While the 
refinement option to move elements or grid reconstruction is attractive since it incurs no extra 
cost, it is the most difficult to implement. In two and three dimensions the difficulties increase. 

As mentioned above, the decay rate a of the spectrum of each element is used as a refine- 
ment criterion to opt between increasing the order of the polynomial, N, and increasing the 
number of elements, K. These are often referred to in the finite element community as p- and 
h-refinement respectively [e.g. 7] and are similar in the spectral element method. 

2.3 Refinement Process 

Once the elements have been marked for refinement and all the limits and tolerances have 
been checked, the refinement process proceeds. There are three types of refinement. 

1) For mesh reconstruction or zero-cost refinement by moving elements and adjusting their 
relative sizes, the elements with the most votes for refinement are shrunk in size according to 
their relative errors with their neighbours. The extra space released by the shrinking of the 
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element is added to the neighbour element, keeping the total space constant. An algorithm is 
implemented to check that the entire domain has been covered and that there are no holes in 
the grid. 

2) For refinement by increasing the order of the polynomial N, the order is simply increased 
by two to N+2. 

3) For refinement by increasing the number of elements K, the element to be refined is 
simply split in two and the polynomial is decreased by two to N-2, with a lower limit imposed 
on the order. The reduction in polynomial order is introduced since the decay rate has shown 
that an Nth order polynomial is not resolving the solution well and is therefore wasteful. 

In two and three dimensions, error estimates and decay rates are obtained for each direction 
in each element and hence the adaptive refinement decisions may be carried out in each direc- 
tion. At all times in the process, the aim is to use the minimum amount of resolution to obtain 
the best solution possible. In the interest of efficiency, there is no “backtracking” to improve 
the solution. Since the problems to be treated are transient this necessitates a tolerance to 
be imposed at all times. The error estimates may be calculated at every step to ensure the 
tolerance is imposed but this can be expensive. Instead, errors are checked periodically with a 
frequency imposed by the user. Since increase in resolution is also expensive, necessitating a 
recalculation of the discrete operator matrix, adaptivity is also done only periodically with a 
frequency imposed by the user. In practice, the no-cost refinement of moving the elements is 
allowed more often than the true refinement procedures of increasing N or K. 


3. Illustrations 


3.1 Burgers’ Equation 


The one-dimensional viscous Burgers equation 


du du d 2 u 

dt + U dx V dx 2 1 


x € [-1; 1], t > 0 


(4a) 


with boundary conditions 

u(-l,i) = tt(M) - 0 (4b) 

and initial conditions 

tz(x,0) = — sxxnrx (4c) 

was chosen to illustrate the capabilities of the adaptive method. For low viscosities, this 
equation admits a nonsingular thin internal boundary layer that must be resolved for spatially 
and temporally accurate numerical solutions to be obtained. An analytical solution is available 
[8,9] and is used to determine accuracy. For u = 0.01/7T the solution develops a sharp gradient 
at the origin at approximately t = l/n. The gradient reaches a maximum at approximately 
t = 0.5. A full study of this equation by several different spectral methods was reported in [10]. 
The conclusion of that study was that spectral methods were not well suited to the calculation 
of thin inner layers. This is especially true if the location of the layers are unknown. The study 
found that accurate results could be obtained but polynomial orders were inordinately high for 
this simple one-dimensional problem. The best spectral method required N=64 for three digit 
accuracy, while the spectral element method obtained four digit accuracy with N = 16 and K— 4. 
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Further, the time step needed to accurately model this transient process was very small due to 
stability restrictions. The present results show that spectral methods are indeed well suited to 
these types of problems provided adaptivity is performed. 

Discretization 

We begin with the temporal discretization. The nonlinear convective term of Burgers’ 
equation is treated explicitly via a third order Adams-Bashforth technique which is stable for 
Courant numbers less than 0.72. The diffusion term is evaluated implicitly via Crank-Nicolson, 
which is unconditionally stable. The time stepping is done in two steps: 


u-u n 1 Jt. d(u n ~ < ') 2 

At ~ 2^°'' dx 

(5a) 



u n +»-* 

= v 4 

At dx 2 

(5b) 

where the a q are the third order Adams-Bashforth coefficients 


23 4 5 

(5c) 

CtQ — —— Qfi = Go = • 

12 3 12 


The spatial discretization follows the standard spectral element formulation [e.g. 4,10,11], 
hence it is only described briefly. The spectral element method breaks up the computational do- 
main into A;=1,K macro-elements of lengths L*, upon each of which the unknown u is expanded 
as a Lagrangian interpolant through the Gauss-Lobatto collocation points: 

«*(»■) = £«?M r ) (6) 

*= 0 


where u 1 - is the 
are given by 


value of the unknown on element k at the collocation point 


h,(r) 


(r* - l)ift(r) 

N(N + l)(r-n)P N (r t ) 


r<E[— 1;1], 


r t . 


The interpolants 


( 7 ) 


Here, Pn is the Nth order Legendre polynomial. The collocation points, r t , are the zeroes of 
the numerator in Equation (7). Each element is mapped to the interval [—1; 1], The discrete 
solution may also be expressed in terms of the spectrum a n as in Equation (1). The equations 
are written in a variational formulation and all integrals are performed by Gauss-Lobatto Leg- 
endre quadrature. 


Results 

The viscous Burgers equation is solved for three cases. For a low viscosity, v — 0 . 01 / 7 T, the 
equation exhibits a sharp gradient as illustrated in Figure la by the spectral element solution 
using N=15 and K— 4. This solution has four digit accuracy in the maximum gradient at the 
origin as shown in Table 1, where it is compared with the analytical and other non-adaptive 
spectral element solutions. This case is used to illustrate the power of the adaptive method 
to refine sharp gradients. For high viscosity, v ~ 1/47T, the diffusion terms dominate and 
the gradient is weak, as shown in Figure lb. This case illustrates how the adaptive method 


5 



behaves in a problem where gradients are weak. The third case involves the solution to Burgers’ 
equation (4) with a new initial condition u = —sin7rx + 0.5 and periodic boundary conditions. 
The solution to this problem is a steepening wave travelling towards x = 1 shown in Figure 2. 
This case illustrates how the adaptive method can resolve non-stationary fronts. 

In each case, the figures show the time evolution of the solution u to Equation (1) and 
the grid used. The grid, represented by symbols, is denoted by the element boundaries (|) 
and the collocation points inside each element (★). Each grid is plotted once when it is newly 
implemented at a vertical axis coordinate corresponding to time — t in Equation (4), except 
where noted. 

1) u = 0.01/tt 

Several properties of the method are illustrated. Each type of refinement will be studied 
separately and the combined adaptive result is shown. 

The first test case, illustrated in Figure 3, shows the merit of refining away from the sharp 
gradients. In Figure 3a, we show a solution where there is sufficient resolution around the sharp 
gradient (N=ll in two elements of size Ax = 0.05 around the origin) but poor resolution away 
from the thin layer (N=3 in four elements of size Ax = .475). Note that the initial condition is 
sufficiently resolved, but as time evolves, errors become significant. The error in the gradient is 
small but the error in time is significant. In Figure 3b, we present the adaptive solution starting 
with the same initial resolution. By adaptively refining the smooth regions by increasing the 
order of the polynomial from N=3 to N = 7 at t=0.1, the errors are reduced significantly. The 
error estimates indicate an error of 10 -3 before refinement and 10 -5 thereafter. 

The second test case, illustrated in Figure 4, shows the merit of no-cost refinement by 
adapting the size and position of elements. In Figure 4a, we show a high resolution solution 
(N=ll, K=4) but with non-optimal equally spaced elements. Note that resolution is adequate 
up until approximately t = 0.3. Beyond, a large oscillation appears indicating poor resolution 
around the sharp gradient. The temporal error is moderate but the error in the maximum 
gradient is significant. In Figure 4b, we present the adaptive solution starting with the same 
initial resolution. By checking every A = 0.5 and adaptively refining by shrinking the elements 
with the greatest errors, we obtain a significantly improved solution. Note that refinement, does 
not start until t = 0.25 as it is not needed until then. At t — 0.45, shortly before the maximum 
gradient develops, the finest grid is chosen with the smallest elements around the sharp gradient 
of size Ax = 0.03125. The tolerance for the error estimates is 10“ 4 . Again the errors in the 
maximum gradient and time at which it occurs are significantly reduced. 

The third test case, illustrated in Figure 5, shows the merit of refining by increasing the 
number of elements. In this case, we choose to illustrate the point with a poor resolution run. 
In Figure 5a, we show a low resolution run with N=3 and K=4. The initial grid consists of 
equally spaced elements. We adaptively refine by moving elements only (no-cost refinement). 
Despite clustering small elements near the sharp gradient, the errors remain large. Note the 
large overshoots. This case shows that refinement by adjusting element sizes is not sufficient 
for improved solutions. In Figure 5b, we show a solution with the same initial grid, where 
refinement has been adaptively performed by increasing the number of elements only. By 
checking the error estimates every A tx = 0.1 and adding elements each time since the errors 
were poor throughout, the adaptive method obtains a vastly improved solution. The final 
grid has K=14 elements and a fixed N=3 per element. The error in the maximum gradient 
is improved from a 66% error to an 11% error. The error in the time at which the maximum 
gradient occurs is also significantly improved. As the errors are still large, it is obvious that 
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none of these refinement methods is optimal. The best refinement procedure combines all three 
as illustrated in the following case. 

The fourth test case is a fully adaptive case where all three refinement methods are used. 
In Figure 6, the time evolution of the solution and of the grid is presented. The element 
demarcations (|) are not shown for clarity. This case was obtained with a frequency of adaptivity 
in element size of A = 0.03, in N of A tpj = 0.075, and in K of A t# = 0.075, an absolute 
tolerance on the error estimators of 5.10“ 5 , a relative tolerance of 5, a maximum on N of 11 
and a maximum on K of 20, a minimum element size of Ax = 0.01 and a minimum time step 
size of 10 -4 . The initial grid is a relatively coarse K=4 equally spaced elements with N=5 in 
each. The final grid has K=16 elements with polynomial orders ranging from N=5 to N=9. 
The low order N=5 elements are generally in the edge areas of the sharp jump and away from 
the jump, while high order polynomials are used in the smooth regions of the jump. Since this 
is a viscous phenomenon, the jump is not a discontinuity but rather a thin layer. Hence it is 
expected that the adaptive algorithm would refine in this way. There are many possibilities 
for an adaptive solution. The choice of adaptivity frequency, absolute and relative tolerances, 
minimum and maximum resolution are left to the user. The results show excellent resolution. 
The maximum gradient is captured with five significant digits and the time at which it occurs 
is accurate to within the smallest time step size used (A£ m j n = 0.0005). Comparing with 
the non-adaptive solutions of Table 1 we see that the solution is better than the corresponding 
N=17 K=4 non-adaptive solution. The greatest expense in a spectral element solution is related 
to the storage (if a direct inversion method is used) and the inversion of the discrete operator 
matrix. The amount of work scales linearly with the number of elements and nonlinearly with 
the polynomial order. Hence the savings afforded by being able to use more elements and lower 
order polynomials is substantial because the bandwidth of the matrix is greatly reduced. In 
this example, the maximum amount of work per iteration presents a savings of 22%. But the 
real savings is associated with the transiently adapting solution. Where the solution has not 
yet developed any fine structures, a coarse grid and a large time step can be used. Including 
these, the savings is of 37% for this example. 

2) v — 1/47T 

For this high viscosity case, the solution exhibits only weak gradients as illustrated in Figure 
lb. We present only one test case, illustrated in Figure 7, to show how the adaptive method 
handles weak gradients. The initial grid has K=4 equal spaced elements with N=5 on each. 
As expected, the adaptive method refines by increasing the order of the polynomial and by 
adjusting element sizes. The error estimators never indicate poor decay rates since this is a 
smooth solution which is ideal for spectral discretization. Hence no elements are ever added in 
the adaptive refinement process. The number of elements remains at K=4 while the polynomial 
order is increased from N=5 to N=7 in the regions of the weak gradient. Note that very little 
resolution is needed in comparison to the sharp gradient case. 

3) Moving Front 

The challenge of the moving front case is to track the steepening wave with sufficient resolu- 
tion, without wasting resolution in areas where the wave has passed through. This case therefore 
depends on a good coarsening algorithm. The results are presented in Figure 8. Again, ele- 
ment demarcations (|) are not shown for clarity. The adaptive method accurately resolves the 
steepening moving front. This case has an initial grid of K=4 equally spaced elements and N = 7 
on each element. The final grid has K = 7 elements with polynomial orders ranging from 7 to 
19. As mentioned above, coarsening is difficult to do without losing accuracy. Coarsening is 


7 



limited to element size and position changes. As a result, the coarsening is somewhat inefficient 
in this case. This is due to the fact that the sharp gradient is adaptively resolved by many 
small elements and that moving these small elements relative to each other is a slow process for 
coarsening. However, the adaptive method successfully tracks the moving front and coarsens 
in the regions behind it once it has moved away. 

The three test cases on Burgers equation show the capabilities of the adaptive method to 
optimally refine solutions in common situations. We now turn to two-dimensional problems 
where adaptivity becomes more complex, but where potential savings increase dramatically. 

3.2 Driven Cavity Flow 

Laminar two-dimensional flow in a cavity is used to illustrate how the adaptive process 
performs in higher dimensions. The full incompressible Navier-Stokes equations 

du 1 

— ^ + -tT-Vu = -Vp+— V 2 w (8a) 

ut K 

V * u = 0, (8b) 

where u is the velocity vector, p is the pressure and R is the Reynolds number, are solved in 
a geometry shown in Figure 9. The flow in the cavity is driven from above with a uniform 
unit velocity. The boundary conditions are all Dirichlet: u=l,t; = 0atj/ = l and no-slip 
boundary conditions on the three walls of the cavity y — 0, x = 0 and x = 1. This flow is a 
popular test case; some references are given in [13]. The discontinuity in velocity at the upper 
corners, where u = 0 from the no-slip condition and u = 1 from the imposed driving flow, can 
present problems in the numerical solution. At these points there is a singularity as the vor- 
ticity becomes infinite. A priori, since the location of the singularity is known, one can resolve 
it adequately in an initial grid. However, the location of singularities is not always known and 
the optimal degree and type of resolution is not known. This example will illustrate how the 
adaptive method refines singularities in a general way. Further, as savings in cpu and storage 
become more important in two and three dimensions, adaptive refinement in two and higher 
dimensions is imperative. The nonconforming spectral element formulation provides the ability 
to locally refine without incurring undue increase in resolution globally. The automation of the 
adaptation greatly enhances the advantages of the nonconforming method as illustrated in this 
example. 

Discretization 

The Navier-Stokes equations (8) are advanced in time using an explicit/implicit fractional 
time-stepping scheme [1,12], which results in a set of separate equations for p and u. The 
nonlinear terms are treated explicitly as a simple inhomogeneity in the elliptic equation set. 
The spatial discretization is a nonconforming spectral element method described in [5]. Here, 
we briefly review the nonconforming formulation. 

The nonconforming discretization allows arbitrary element matchups in the computational 
grid, an example of which is shown in Figure 10b. While C° continuity is lost due to the mis- 
match of approximating polynomials on either side of an elemental interface, the consistency 
error is held as small as the approximation errors, thereby preserving the convergence proper- 
ties of the spectral element discretization. The implementation consists of introducing a new 
structure known as mortars 7 P , which are defined as the intersection of adjacent element edges. 
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Upon this structure, </> is defined as the mortar function, which is a polynomial of degree N 
in the local one-dimensional mortar variable s. The approximation space X h consists of the 
functions v in C 2 that are tensor products of polynomials of degree N in each direction of each 
element A;, such that the two following conditions are satisfied: 

1- the vertex condition: at each vertex q of each element Ar, v\ k (q) = (j>(q ). 

2- the C 2 condition: over each elemental edge, 6 Pn - 2 f e dge( v ~ 4 > )' l P^ s = 0- 

For a conforming approximation this definition of Xh reduces to the standard spectral element 
approximation space. The vertex condition ensures exact continuity at cross points while the 
C 2 condition represents a C 2 minimization of the jump in functions at internal boundaries. The 
combination of these two conditions ensures the optimality of the discretization, as explained 
and illustrated in [5]. 

Results 

The driven cavity flow is solved adaptively starting with a coarse N=4 K— 4 equal-sized 
element grid shown in Figure 10a. In each figure, the grid is shown by solid element boundary 
lines. The collocation points are not shown. The Reynolds number for this test case is #=100. 
The streamlines are shown for purposes of comparison. In this test case, refinement is limited to 
moving elements, splitting elements and increasing the polynomial order globally. 1 he option 
to increase order locally is not used. For this steady-state problem, adaptivity in time is not 
crucial. Figure 10 presents the intermediate adaptation steps. The first adaptation (Figure 
10b) subdivides the upper corner elements into four as both x and y direction decay rates are 
below 1. The number of elements grows from K=4 to K = 10. The second adaptation (figure 
10c) shrinks the corner elements at the expense of their neighbours. This is a no-cost refine- 
ment step which decreases errors only slightly. The third adaptation (not shown) increases the 
order of the polynomial globally to N— 6. The fourth adaptation (Figure lOd) subdivides the 
corner elements in the y direction. The number of elements grows from K — 10 to K = 12, while 
the polynomial order remains N=6 from the previous step. Because the singularity is difficult 
to resolve with smooth polynomials, the error estimator always indicates poor decay rates in 
the elements containing the singularities. Hence, refinement by adding or spplitting elements 
is dominant in this case. In two dimensions, as in one, the error estimators perform very well. 
The decay rate is measured in both x and y directions, providing adaptivity criteria for both 
directions. In adapting from the mesh in Figure 10c to that in Figure lOd for example, the 
decay rate is below 1 only in the y direction for the corner elements; hence adaptivity proceeds 
by splitting the corner elements vertically only. The error estimates are reduced from 10 1 
to 10 -4 in most of the domain. In the corner elements containing the singularities, the error 
estimates are reduced to 10“ 2 . These relatively poor errors are expected as the singularity is 
hard to resolve with smooth functions. Further refinement could reduce this error. However, 
the main advantage of the adaptive method is demonstrated: the error due to the singularity 
is localized to a very small element and contamination of the rest of the domain is minimized. 


4. Conclusions 

The test cases presented clearly show that the adaptive formulation of the spectral element 
method increases the flexibility and capabilities of the method. Sharp gradients, singularities 
and regions of poor resolution can be resolved optimally. The method combines three types of 
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refinement using criteria based on error estimators which indicate the quality of the resolution 
on each element. The savings in storage and cpu time appear to be considerable. In two and 
three dimensions the savings increase. The increase in efficiency should broaden the range 
of problems that can be investigated by high order spectral type methods. Generally, the 
savings is afforded by using a larger number of elements, with very high order only in elements 
where it is really needed and effective. These are regions of smooth solutions. In regions 
of discontinuities or sharp structures many lower order elements are used. A larger number 
of elements can also increase the efficiency of coarse-grained parallel implementation of the 
method [14], since a larger number of elements can be distributed among many processors. 
The overhead associated with adaptivity, however, can be substantial. The error estimators are 
relatively cheap to calculate but the inversion of the matrix that is needed with each new grid 
is expensive. This expense is minimized by adapting only periodically. While the refinement 
capabilities were successfully demonstrated here, the coarsening algorithm remains somewhat 
inefficient. A more clever algorithm will be the subject of future work. 

This study of the practical aspects of adaptive meshes serves as a learning step in the on- 
going development of a fully automatic adaptive mesh procedure for spectral element methods 
in two and three dimensions. The purpose of this work is to make high order spectral methods 
more accessible to engineers for the simulation of complex physical phenomena by automating 
the grid generation and refinement processes and by reducing resolution requirements signifi- 
cantly. 
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Figure 1: Time evolution of the solution to the viscous Burgers equation with homogeneous 
bounday conditions and sine wave initial condition a) low viscosity v = 0.01/7T sharp gradient 
case b) high viscosity v — 1/4tt weak gradient case 
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Figure 3: Effect of refinement by increasing polynomial order in smooth regions of the solution to 
the viscous Burgers equation for low viscosities: a) fixed resolution solution b)increased resolution 
solution 



Figure 4: Effect of no-cost refinement by adapting element size and position for the solution to 
viscous Burgers equation for low viscosities a) fixed equally-spaced element solution b) adaptive 
solution 





Figure 5: Effect of refinement by adapting increasing the number of elements for the solution to 
viscous Burgers equation for low viscosities a) fixed low order but adaptive element size solution 
b) adaptive (number of elements) solution 
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Figure 6: Fully adaptive solution to viscous Burgers equation for low viscosities: final grid has 
K=16 and N varying from 5 around the leading and trailing edges of the sharp jump to N=9 
within the jump 



Figure 7: Fully adaptive solution to viscous Burgers equation for high viscosity: v — 1/47T 
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Figure 8: Fully adaptive solution to viscous Burgers equation for moving front case: both refine- 
ment and coarsening are used 


U=1 y=Q 

11=0 v=0 

u=0 v=0 

Figure 9: Geometry for two-dimensional Navier-Stokes fluid flow in a driven cavity 
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